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Lattice and Implications for Bose-Hubbard Model 
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We show that one can properly take into account of the interaction effects and construct a set of or- 
thonormal Wannier functions for a Bose-Einstein condensate in an optical lattice. These interaction- 
dependent Wannier functions are used to compute the tunneling rate J and the on-site repulsion 
U in the Bose-Hubbard model. Both parameters are found to be substantially different from ones 
calculated with the single-particle Wannier functions. Our numerical results of U are found in good 
agreement with the measured on-site energy in a recent experiment [Campbell et al. Science 314, 
281 (2006)]. 
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The system of a Bose-Einstein condensate (BEC) in 
an optical lattice has been one of the most exciting and 
studied systems in recent years P, 0, S S|- The sys- 
tem has enabled experimentalists to observe for the first 
time many interesting phenomena predicted a long time 
ago in condensed matter physics. The most spectacular 
example is the observation of the quantum phase transi- 
tion from a superfluid to a Mott insulator in such a BEC 
system 0,(1 [11,0. 

To understand this periodic BEC system, one often 
uses the Wannier functions and reduces the system to 
the famed Bose-Hubbard model (BHM) [TO, M, H El • 
The Hamiltonian of the BHM is given by 



(1) 



where hi and a] are respectively the bosonic annihilation 
and creation operators at the ith lattice site and hi = 
a\a,i. The angle brackets above indicate the summation 
over all nearest neighboring pairs. The tunneling rate J 
and the on-site repulsion U can be computed with the 
Wannier function w(r) [ll|, E2, EH ■ For atoms of mass m 



and s-wave scattering length a s , they are given by 11 1 
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where and Vj are the coordinates of a pair of nearest 
neighboring sites. The potential for the optical lattice has 
the form Vi att = VoE R Y^f=i sin 2 (qBXj) with q B = n/d 



being the laser wave vector, d the lattice period and Vo 
the laser intensity in units of the recoild energy En = 
h 2 q 2 B /2m. D = 1, 2, 3 is the dimensionality of the lattice. 
As is well known [13|, when the filling factor, namely 
the average number (hi) of bosons at one site, is fixed, 
the physics of the BHM is completely determined by its 
tunneling rate J and on-site repulsion U. Hence it is 
crucial to have accurate values of J and U for a good 
description of the BEC system with the BHM. 

So far, these two parameters are usually computed 
with the single-particle Wannier function [ll|, |l2j, Il3 |. 
Such treatment is good only in the low-filling regime [3|. 
For higher fillings, due to stronger inter-atomic interac- 
tions, the shape of the Wannier function is expected to 
distinguish significantly from that of the single-particle 
Wannier function. As a result, J and U are interaction- 
dependent. In particular, J is more sensitive since it 
depends on the tails of the Wannier function as seen in 
Eq. ([2])- This view is echoed in literature. For example, 
Bloch et al. pointed out [3j, "For intermediate fillings, 
the Wannier functions entering both the effective hopping 
matrix element J and on-site repulsion U have to be ad- 
justed to account for the mean-field interaction". In a 
recent experiment by Campbell et al. the one-site 

energy with the filling factors increasing from one to five 
was observed to have a 27% decrease. 

To our best knowledge, there has been only one sys- 
tematic theoretical attack on this important problem. 
This was carried out by Li et al. in Ref. [151], where 
the authors constructed a set of orthonormal interaction- 
dependent Wannier functions with Kohn's variational ap- 
proach [lij ]. This variational method has one intrinsic 
shortcoming: the chosen trial Wannier function may be 
quite different from the true Wannier function. In Ref. 
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|15j ]. Li et al. used Gaussian functions as the trial func- 
tions. However, as pointed out in Ref. 0, [l3|, the Gaus- 
sian approximation is not good for calculating the tunnel- 
ing parameter J. This means that even the best Gaussian 
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type Wannier functions obtained by variation in Ref. 
are not good enough. 

In this Letter we show that a set of Wannier func- 
tions can be constructed from the Bloch states of the 
Gross-Pitaevskii equation (GPE) for the periodic BEC 
system. These Wannier functions are proved to be or- 
thonormal and are therefore suitable for the use of redu- 
ing the BEC system to the BHM. Moreover, these Wan- 
nier functions are interaction-dependent by construction; 
we call them nonlinear Wannier functions to distinguish 
from the single-particle Wannier functions. That the in- 
teraction effects are properly taken into account in these 
nonlinear Wannier functions roots in the fact that the 
Bloch states used to construct them minimizes the sys- 
tem energy under the mean-field approximation. With 
these nonlinear Wannier functions, we find that both the 
tunneling rate J and the on-site repulsion U are substan- 
tially affected by the mean-field interaction. For simplic- 
ity, the construction and properites of the nonlinear Wan- 
nier functions are illustrated in detail for one-dimensional 
optical lattice while the results for the BHM parameters 
J and U are presented for all dimensionality. 

To define the nonlinear Wannier functions, we first in- 
troduce the mean- field Bloch states ipk{ r ) [3 fl9l |2Cl| . 
which satisfy the following time-independent GPE [2l| 
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where g = Airh 2 a s /m and no is the average BEC den- 
sity. The band index is omitted as we focus on the low- 
est Bloch band. The nonlinear Wanniner functions are 
constructed from these Bloch states as follows 
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where the integration is over the first Brillouin zone and 
SI | is its volume. This is exactly the same way how 
the single-particle Wannier function is constructed from 
Because of the nonlinear term 
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Bloch states 

in Eq. (g]), one may doubt whether the Wannier func- 
tions constructed in such a way are orthonormal to each 
other and thus have any use. This doubt can be cast 
aside immediately by observing that the Bloch states 
defined in Eq. (J4| are orthonormal to each other, i.e., 



/ V£ (r)^ k (r)dr 



J k'k' 



despite the nonlinear term in 



Eq. (|4]). With this, one can prove the orthonormality of 
the nonlinear Wannier functions 



J w*(r — Ti)w(r — Tj)dr 
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where the integration is over the entire space. There was 
a concern in Ref. ,15] that the definition in Eq. © 



would fail because of the existence of a loop structure 
in the mean-field Bloch band /x(k) [H, O, 0]. This 



concern is not justified because the BHM is a single-band 
approximation and it is a good description of the BEC 
system only when the band gap of the lattice potential 
is much larger than the interatomic interaction 11| . The 
loop structure in the Bloch band appears only for shallow 
optical lattices, where the BHM does not apply. 

We now use the ID case to illustrate some key points 
in the numerical computation of the nonlinear Wannier 
functions. The complete set of the Bloch functions ipk{%) 
in the lowest Bloch band can be obtained by solving Eq. 
dU with the same method in Ref. [251 ] . Since an arbitrary 
phase can be added to each ipk(x), to obtain a proper 
Wannier function via Eq.©, the numerical method must 
be designed to make sure that the resulted ?pk (x) is ana- 
lytic in k. Usually, one also wants the Wannier function 
to be real and symmetric (or antisymmetric) . To achieve 
this, one should pay attention to the values of ipo (0) and 
V-Vd(O) (i) If both ipo(0) and ip w /d(0) are nonzero, 

the phase of the Bloch function must be chosen such that 
■0fc(O) is real, (it) If both ^o(O) and ipn/di^) vanish, the 
phase must be chosen such that V'fe(O) is purely imagi- 
nary, (in) If only one of ^o(O) and V^/dCO) is zero, one 
can shift the origin in the x space by half of the lattice 
constant. With the new origin, one is then back to either 
case (i) or (ii). One can prove that (1) if ^o(O) ^ 0, 
w(— x) = w(x) and w*(x) = w(x), that is, the Wannier 
function w(x) is symmetric about x = and real; (2) if 
"0o(O) = 0, w(x) is antisymmetric about x = and real. 
Following Kohn's strategy [III], one can also prove that 
no other choices of phases in ipk(x) can lead to nonlin- 
ear Wannier functions that are both real and symmetric 
(antisymmetric) about x = 0. 
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FIG. 1: (color online) ID nonlinear Wannier functions w(x) 
for different mean-field interactions gno. The insets show the 
decay of the local maximums of the ID nonlinear Wannier 
function. 

Our numerical results of the ID nonlinear Wannier 
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functions are plotted in a semilog fashion in Fig. [TJ As 
clearly shown, the decay of the the Wannier function re- 
mains exponential despite the nonlinearity. The effect 
of the nonlinearity (or interaction) is to make the de- 
cay slower and thus the Wannier function less localized. 
The single-particle Wannier function has been proven to 
decay exponentially by Kohn[22|. It was pointed out 
later [26| that this exponential decay is related to a well- 
known mathematical result that connects the behavior of 
a function near a branch point to the asymptotic decay 
of its Fourier transform [U [Hj]. Suppose that fk(x) 
is a periodic function fk(x) — fk+2n/d( x ) an d has a 
leading behavior at the branch point fco = tt / d + ih as 

fk( x ) = fo( x ) + 7 [i{k — ko)] 13 . Its Fourier transform is 
an exponential decay function, F(x) — J dkfk(x)e~ tkx — 
27sinvr(l + 7)r(l+7)2;"( 1+7 )e-' l:c . This means that the 
nonlinear Wannier functions, shown to decay exponen- 
tially in Fig. [TJ may also have these analytical properties. 
The rigorous proof for this, however, is left for the future 
investigation. Also note that due to the power-law pref- 
actor [28| , the decay in Fig. [1] is not strictly exponential 
as indicated by the slight curving of the envelope of the 
peaks. 

Nonlinear Wannier functions can be computed simi- 
larly for the 2D and 3D optical lattices. However, the 
computation time can become enormous in particular for 
the 3D case. Since our ultimate goal is to compute the 
two basic parameters J and U of the BHM from Eqs. 
([2j) and respectively, we can reduce the computa- 
tion time significantly by not calculating out the Wan- 
nier function explicitly. For J, one can combine Eq.(|2j) 
and Eq.© to express J in terms of the Bloch functions 
and then use this expression to compute J efficiently. 
For U, one can reduce the computing time by utilizing a 
well-known fact that U is proportional to the difference 
between /j,(k = 0) and the system's mean-field ground 
state energy. 

Our numerical results for J and U are shown in Figs [2] 
and[3]for different values of Vo and gno- As clearly shown 
in the figures, the mean-field interactions gno have pro- 
nounced effects on both the tunneling rate J and on-site 
repulsion U. For a fixed lattice strength Vo, with the 
increase of gno, J increases while U decreases dramati- 
cally. This is expected as the Wannier function becomes 
less localized as gno increases. Moreover, it is apparent 
from the figures, in higher dimensions, the interaction ef- 
fects are much stronger. This highlights the need to take 
into account the interaction effect into J and U since 
most of the experiments are carried out with the 3D op- 
tical lattices. Since the optical lattices have the form 
Viatt = VqEr J2f=i sm2 (iBXj), the different directions 
in our system are decoupled in the linear case gno = 0. 
The dependence of the interaction effects on dimension- 
ality shows that the interaction can strongly couple the 
motions along different directions. 

Let us take a closer look at a BEC in a 3D optical lat- 
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FIG. 2: (color online) Tunneling rate J via the mean-field in- 
teraction gno for optical lattices with different strength Vo. 
ID, 2D, and 3D correspond to the dimensionality of the op- 
tical lattice. The Tunneling rate J, mean-field interaction gno 
and lattice depth Vo are all in units of Er. 
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FIG. 3: (color online) On-site repulsion U via the mean-field 
interaction gno for optical lattices with different strengths 
Vo. ID, 2D, and 3D correspond to the dimensionality of 
the optical lattice. Repulsive on-site interaction U is in units 
of Enas/d; the mean-field interaction gno and lattice depth 
Vo are in units of Er. 



tice with Vq = WEr, At a low filling with gno = O.OIEr, 
we have J/E R = 0.019 and Ud/E R a s = 22.4 (see Figs. 
[2](b) and[3](b)), which is consistent with the results cal- 
culated with the single-particle Wannier function in Ref. 
[2!j[3(3]. However, at a higher filling with gno — 2.0Er, 
J is approximately three times larger than the value cal- 
culated with the single-particle Wannier function. Mean- 
while, the on-site energy U d/E R a s is significantly re- 
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duced to 6.6. This strong dependence of J and U on 
the grig justifies the necessity of introducing the nonlin- 
ear Wannier functions. In this work, the parameters gn^ 
in Figs. [2] and [3] relates to the filling factor (rii), which is 
often used in the literature @, as n = (rn) /d 3 . 

In a recent experiment by Campbell et al. 14], the 
on-site energy U was measured for different filling fac- 
tors with the two-photon Bragg spectroscopy [3l| . They 
found that U = 22Hz for the (n,) = 5 shell at V = 35E R , 
a decrease of 27% from U = 30Hz for (rij) = 1 ~ 2. 
Our numerical results are U = 28.4Hz for (n,) = 5 and 
U = 33.1 ~ 31.8Hz for (n,) = 1 ~ 2, in good agreement 
with the experiment. This shows that even though the 
nonlinear Wannier function is a mean-field concept, it 
somehow still captures much of the interaction effect in 
the Mott insulator regime. A possible method for the ex- 
perimental study on the tunneling rate J is to study the 
interference pattern produced by an expanding atomic 
cloud 0. 

We emphasize here that our calculation of the two 
basic parameters of the BHM has been done with the 
mean-field theory. Further improvement of the theoret- 
ical framework is also needed to include the effects of 
quantum fluctuations [32I ]. 

To conclude, we have demonstrated a way to construct 
a set of orthonormal Wannier functions by properly incor- 
porating interaction effect for a BEC in an optical lattice. 
Although these nonlinear Wannier functions are less lo- 
calized than the single-particle Wannier functions due to 
the repulsive interaction, they retain many of the analyt- 
ical properties of the single-particle Wannier functions. 
For example, our numerical results show that they decay 
exponentially. We have used these Wannier functions to 
compute the tunneling rate J and the on-site interaction 
U in the Bose-Hubbard model. The computed U are 
found in good agreement with experimental results. 
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